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We address the thermodynamics, density profiles and superfluid density of trapped fermions undergoing BCS- 
BEC crossover. We consider the case of zero and finite population imbalance. Our approach represents a fully 
consistent treatment of "pseudogap effects". These effects reflect the distinction between the pair formation 
temperature T* and the pair condensation temperature Tc- As a natural corollary, this temperature difference 
must be accommodated by modifying the fermionic excitation spectrum iJk to reflect the fact that fermions 
are paired at and above Tc- It is precisely this natural corollary which has been omitted from all other many 
body approaches in the literature. At a formal level, we show how enforcing this corollary implies that pairing 
fluctuation or self energy contributions enter into both the gap and the number equations; this is necessary in 
order to be consistent with a generalized Ward identity. At a less formal level, we demonstrate that we obtain 
physical results for the superfluid density ris (T) at all temperatures. In contrast, previous work in the literature 
has led to non-monotonic, or multi- valued or discontinuous behavior for ria (T). Because it reflects the essence 
of the superfluid state, we view the superfluid density as a critical measure of the physicality of a given crossover 
theory. In a similarly unique fashion, we emphasize that in order to properly address thermodynamic properties 
of a trapped Fermi gas, a necessary first step is to demonstrate that the particle density profiles are consistent with 
experiment. Without a careful incorporation of the distinction between the pairing gap and the order parameter, 
the density profiles tend to exhibit sharp features at the condensate edge, which are never seen experimentally 
in the crossover regime. The lack of demonstrable consistency between theoretical and experimental density 
profiles, along with problematic behavior found for the superfluid density, casts doubt on previous claims in the 
literature concerning quantitative agreement between thermodynamical calculations and experiment. 



PACS numbers: 03.75.Hh, 03.75.Ss, 74.20.-z 



I. INTRODUCTION 



There has been a resurgence of interest in studies of the 
crossover between the usual BCS form of fermionic super- 
fluidity and that associated with Bose Einstein condensation 
(EEC). This is due, in part, to the widespread pseudogap phe- 
nomena which have been observed in high temperature su- 
perconductors in conjunction with the small pair size. The 
latter, in particular, was argued by Leggett to be a rationale' 
for treating the cuprates as mid-way between BCS and BEC. 
Others have argued^i^ that the cuprate pseudogap can be un- 
derstood as arising from pre-formed pairs which form due to 
the stronger-than-BCS attraction. Additional reasons for the 
interest in BCS-BEC crossover stem from the precise realiza- 
tion of this scenario in ultracold trapped Fermi gases, where 
the attractive interaction can be continuously tuned from weak 
to strong via a Feshbach resonance in the presence of a mag- 
netic field. A final rationale for interest in this problem stems 
from the fact that BCS theory is the prototype for successful 
theories in condensed matter physics; and we now have come 
to realize that this is a very special case of a much more gen- 
eral class of superfluidity. 

BCS-BEC crossover theory is based on the observation^-^ 
that the usual BCS ground state wave function 'fo — nk(uk + 
VkC^ |clj. |)|0) (where cj^ ^ and c^. ^ are the creation and an- 
nihilation operators for fermions of momentum k and spin 
a =t, i) is far more general than was initially appreciated. If 
one tunes the attractive interaction from weak to strong, along 
with a self consistent determination of the variational param- 
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eters and Uk the chemical potential passes from positive 
to negative and the system crosses continuously from BCS to 
BEC. The vast majority (with the possible exception of the 
high Tc cuprates) of metallic superconductors are associated 
with weak attraction and large pair size. Thus, this more gen- 
eralized form of BCS theory was never fully characterized or 
exploited until recently. There are a number of different ren- 
ditions of BCS-BEC crossover theory. Each rendition can be 
represented by a selected class of many-body Feynman dia- 
grams, often further simplified by various essential or non- 
essential approximations. There is no controlled small param- 
eter and thus the selection process is based on highly variable 
criteria. For the most part the success or failure of a particular 
rendition is evaluated by comparing one or a set of numbers 
with experiment. 

It is the goal of the present paper to discuss a criteria set for 
evaluating BCS-BEC crossover theories which captures the 
crucial physics, rather than the detailed numerics. We apply 
these criteria successfully to one particular version of BCS- 
BEC crossover theory which builds on the above ground state. 
In this context we address a wide range of physical phenom- 
ena. These include density profiles, thermodynamical prop- 
erties and superfluid density with application to polarized as 
well as unpolarized gases. It is our philosophy that appropri- 
ate tests of the theory should relate to how qualitatively sound 
it is before assessing it in quantitative detail. Detailed quan- 
titative tests are essential but if the qualitative physics is not 
satisfactory, quantitative comparisons cannot be meaningful. 

Four important and inter-related physical properties are em- 
phasized here, (i) There must be a consistent treatment of 
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"pseudogap" effects. As a consequence of the fact that the 
pairing onset temperature T* is different"'^ from the conden- 
sation temperature T^, the fermionic spectrum, E\s, must nec- 
essarily reflect the formation of these pairs. To accommodate 
the pseudogap, i?k must be modified from the strict BCS form 
which has a vanishing excitation gap at and above Tc. Every- 
where in the Hterature an unphysical form for is assumed 
except in our own work and briefly in Ref. 8. (ii) The theory 
must yield a consistent description of the superfluid density 
ng{T) from zero to T^.. The quantity ns{T) should be single 
valued, monotonic,^ and disappear at the same one com- 
putes from the normal state instability. Importantly, Ug (T) is 
at the heart of a proper description of the superfluid phase, 
(iii) The behavior of the density profiles, which are the basis 
for computing thermodynamic al properties of trapped gases, 
must be compatible with experimental measurements. Near 
and at unitarity, and in the absence of population imbalance, 
they are relatively smooth and featureless, unlike a true BEC 
where there is clear bimodality. This can present a challenge 
for theories which do not accommodate pseudogap effects and 
which then deduce sharp features at the condensate edge, (iv) 
The thermodynamical potential fl should be variationally con- 
sistent with the gap and number equations. It should satisfy 
appropriate Maxwell relations and at unitarity be compatible 
with the constraint"^ " relating the pressure p to the energy. 

There has been widespread discussion about the role of 
collective modes in the thermodynamics of fermionic super- 
fluids. And this has become, in some instances, a basis for 
additional evaluation criteria of a given BCS-BEC crossover 
theory. Because the Fermi gases represent neutral superflu- 
ids with low lying collective modes, one might have expected 
these modes to be more important than in charged supercon- 
ductors. Nevertheless, the BCS wave function and its associ- 
ated finite temperature behavior is well known to work equally 
well for charged superconductors and neutral superfluids such 
as helium-3. In strict BCS theory thermodynamical properties 
are governed only by fermionic excitations. This applies as 
well to the superfluid density (in the transverse gauge). Col- 
lective modes are important in strict BCS theory primarily to 
establish that (T) is properly gauge invariant. 

One can arguei? that collective modes should enter ther- 
modynamics as the pairing attraction becomes progressively 
stronger. The role of these modes at unitarity is currently un- 
resolved. In the Bogoliubov description of a true Bose super- 
fluid there is a coupling between the pair excitations and the 
collective modes, which results from inter-boson interactions. 
Thus it is reasonable to expect that the collective modes are 
important for thermodynamical properties in the BEC regime. 
At the level of the simple mean BCS-Leggett wave function 
we find that, just as in strict BCS theory, the collective modes 
do not couple to the pair excitations; this leads to a /2M* 
form of the pair dispersion. The low-lying collective mode 
dispersion is^^ of course, linear in q. All inter-boson effects 
are treated in a mean field sense and enter to renormalize the 
effective pair mass M* . To arrive at a theory more closely 
analogous to Bogoliubov theory, one needs to add additional 
terms to the ground state wave function- consisting of four 



and six creation operators in the deep BEC. The complexity 
becomes even greater in the unitary regime, and there is, in 
our opinion, no clear indication one way or the other on how 
the pair excitations and collective modes couple. 

Our rationale for considering the simplest ground state 
wave function (which minimizes this coupling) is as follows. 
It is the basis for zero temperature Bogoliubov-de Gennes 
(BdG) approaches which have been widely applied to the 
crossover problem. It is the basis for a T = Gross-Pitaevskii 
description in the far BEC regime.'^ It is the basis for the bulk 
of the work on population imbalanced gases. At unitarity the 
universality relation'" " between pressure and energy holds - 
separately for the fermionic contribution (which is of the usual 
BCS form with an excitation gap distinct from the order pa- 
rameter) and for the bosonic term, due to the q^ form of the 
pair dispersion. Finally, this wave function is simple and ac- 
cessible. Thus, it seems reasonable to begin by addressing the 
finite T physics which is associated with this ground state, in 
a systematic way. 

The remainder of this paper presents first the theoretical 
framework for the principal self consistent equations describ- 
ing the total excitation gap, the order parameter, and the num- 
ber equation or fermionic chemical potential. The conse- 
quences for thermodynamics, density profiles and the super- 
fluid density are then presented in separate sections, along 
with numerically obtained results for each property. We 
discuss these properties at the qualitative as well as semi- 
quantitative level, in the context of comparison with experi- 
ment. In the conclusions section, we present a summary of 
the strengths and weaknesses of the present scheme. 



II. THEORETICAL BACKGROUND 

A. Early Relevant History of BCS-BEC Crossover 

While the subject began with the seminal T = work by 
Eagles^ and Leggett,'' a discussion of superfluidity beyond 
the ground state was first introduced into the literature by 
Nozieres and Schmitt-Rink.'^ Randeria and co-workers refor- 
mulated this approach'^ and moreover, raised the interesting 
possibility that crossover physics might be relevant to high 
temperature superconductors^ Subsequently other workers 
have applied this picture to the high Tc cuprates'^''^''** '"^ and 
ultracold Fermi gases^*^'^' as well as formulated alternative 
schemes22i2,? for addressing T 7^ 0. 

The recognition that one should distinguish the pair for- 
mation temperature T* from the condensation temperature 
Tc was cruciali^ii Credit goes to those who noted that pseu- 
dogap effects would appear in the BCS-BEC crossover sce- 
nario of high temperature superconductors, notably first in the 
spin channel.^"* Shortly thereafter, it was recognized that these 
important pseudogap phenomena also pertain to the charge 
channel i'^-^^-^^ And finally, we make note of those papers 
where the concept of pseudogap effects was introduced into 
studies of the ultra-cold gasesj^i2i22 
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B. Pair Fluctuation Approaches to Crossover 

In this section we discuss the present scheme for BCS-BEC 
crossover, as well as compare it with alternative approaches 
including that of Nozieres and Schmitt-Rink. The Hamil- 
tonian for BCS-BEC crossover can be described by a one- 
channel model. In this paper, we address primarily a short 
range s-wave pairing interaction, which is often simplified as 
a contact potential US{x — x'), where U < 0. This Hamil- 
tonian has been known to provide a good description for the 
crossover in atomic Fermi gases which have very wide Fes- 
hbach resonances, such as and ^Li. The details are pre- 
sented elsewhere.^ 

We begin with a discussion of T-matrix based theories. 
Within a T-matrix scheme one considers the coupled equa- 
tions between the particles (with propagator G) and the pairs 
[which can be represented by the T-matrix t{Q)] and drops 
all higher order terms. Without taking higher order Green's 
functions into account, the pairs interact indirectly via the 
fermions, in an averaged or mean field sense. The propaga- 
tor for the non-condensed pairs is given by 

t;giQ) = u-' +xiQ), (1) 

where U is the attractive coupling constant in the Hamilto- 
nian and x is the pair susceptibility. The function x{Q) is the 
most fundamental quantity in T-matrix approaches. It is given 
by the product of dressed and bare Green's functions in var- 
ious combinations. One could, in principle, have considered 
two bare Green's functions or two fully dressed Green's func- 
tions. Here, we follow the work of Ref. 28. These authors 
systematically studied the equations of motion for the Green's 
functions associated with the usual many body Hamiltonian 
for superfluidity and deduced that the only satisfactory trun- 
cation procedure for these equations involves a T-matrix with 
one dressed and one bare Green's function. The presence of 
the bare Green's function in the T-matrix and self-energy is 
a general, inevitable consequence of an equations of motion 
procedure. 

In this approach, the pair susceptibility is then 

x{Q)^Y.^o{Q-K)G{K), (2) 

K 

where Q = {iili,q), and G and Go are the full and bare 
Green's functions respectively. Here G'f^^{K) ~ iujn — ^k, 
Ck = — Ck — h'^k'^ /2m is the kinetic energy of fermions, 
and fi is the fermionic chemical potential. Throughout this pa- 
per, we take h ^ 1, ks = 1, and use the four- vector notation 
K EE (iw„,k), Q = (ir2;,q), J2k = ^E„Ek' etc, where 
u!n — {2n + l)7rT and 51; = 21ttT are the standard odd and 
even Matsubara frequencies^'^ (where n and / are integers). 
The one-particle Green's function is 

G-\K) - ILUr, - - ^{K) , (3) 

where 

E(i^)=^i(Q)Go(Q-/0, (4) 
Q 



More generally, either Go or the fully dressed G is introduced 
into ), according to the chosen T-matrix scheme. Finally, 
in terms of Green's functions, we readily arrive at the number 
equation: n = J2K,a Ga{K). 

Because of interest from high temperature superconductiv- 
ity, alternate schemes, which involve only dressed Green's 
functions have been rather widely studied. In one alternative, 
one constructs a thermodynamical potential based on a cho- 
sen self-energy. Here there is some similarity to that T-matrix 
scheme which involves G only. One variant of this "conserv- 
ing approximation" is known as the fluctuation exchange ap- 
proximation (FLEX) which has been primarily applied to the 
normal state. In addition to the particle-particle ladder dia- 
grams which are crucial to superfluidity it also includes less 
critical diagrams in the particle-hole channel; the latter can 
be viewed as introducing spin correlation effects. Since it in- 
volves only dressed Green's functions, one evident advantage 
of this approach is that it is "I>-derivable^' or conserving. This 
implies that because it is based on an analytical expression 
for the thermodynamical potential, thermodynamical quanti- 
ties obtained by derivatives of the free energy are identical to 
those computed directly from the single particle Green's func- 
tion. 

For a variety of reasons this FLEX scheme, as applied to 
superfluids and superconductors, has been found to be prob- 
lematic. The earliest critique of the GG, T-matrix scheme is 
in Ref. 28. The authors noted that using two dressed Green's 
functions "could be rejected by means of a variational princi- 
ple". They also observed that there would be an unphysical 
consequence: a low T specific heat which contained a con- 
tribution proportional to T^. In a related fashion it appears 
that the FLEX or GG, T-matrix scheme is not demonstrably 
consistent with the Hamiltonian-based equations of motion. 
There also is concern that considering only dressed fermion 
propagators, G, may lead to double counting of Feynman dia- 
grams. Vilk et aP^ noted that the FLEX scheme will not pro- 
duce a proper pseudogap, due to the "inconsistent treatment 
of vertex corrections in the expression for the self energy." 

By dropping the non-dominant particle-hole diagrams, oth- 
ers have found a more analytically tractable scheme^-^. How- 
ever, this scheme fails to yield back BCS-like spectral prop- 
erties which would be anticipated above Tc in a BCS-BEC 
crossover scenario. Among the unusual features found is a 
four excitation branch structure, ■^"^■^'^ not compatible with the 
expected pseudogap description, which should reflect precur- 
sor superconductivity effects in the normal state. In this pseu- 
dogap picture,^ there would be two peaks in the spectral func- 
tion, rather than four More recently, the authors of Ref. [1^ 
applied a related conserving approximation below Tc- They 
did not consider particle-hole diagrams, but included in the 
particle-particle channel a "twisted" ladder diagram. These 
authors found that there was a discontinuity in the transition 
temperature calculated relative to that computed-'^ above Tc- 
They, then, inferred that at unitarity there is a first order phase 
transition, which has not been experimentally observed. 

In the NSR scheme, which is, perhaps, the most widely 
applied of all pair fluctuation theories, one uses two bare 
Green's functions in xiQ) for the normal state. Within this 
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NSR approach, the results are generally extended below Tc by 
introducing^*^ into x{Q) the diagonal and off-diagonal forms 
of the Nambu-Gor'kov Greens functions. At the o utset, the 
fermionic excitation spectrum i?k = v^k+^sc involves 
only the superfluid order parameter. Age, so that the fermions 
are treated as gapless at and above Tc, despite the fact that 
there is an expected "pseudogap" associated with pairing on- 
set temperature T* . The original authors'^ suggested that pair 
fluctuations should enter into the number equation, but ap- 
proximated their form based on only the leading contribu- 
tion in the Dyson series. This approximate form was intro- 
duced via contribution to the thermodynamical potential Q.. A 
more systematic approach, which is based on a full Dyson re- 
summation leads to a form equivalent to Eq. (|4]i, with a bare 
Xo{Q) = TliK ^o{K)Go{Q — K), as was first pointed out 
in R ef. |39i . This more complete scheme was implemented in 
Ref-di 

Another important aspect of the NSR scheme should be 
noted. Because the pairing fluctuation contributions do not 
enter into the gap equation, the gap equation cannot be de- 
termined from a variational condition on the thermodynamic 
potential. In this regard, a rather different alternative to the 
approximated number equation of Ref. was recently intro- 
duced in Ref. l40pll These authors argued one should com- 
pensate for the fact that dfl/dAgc 7^ by adding a new term 
(deriving from this discrepancy) to the number equation. We 
view this latter alternative as even more problematic since it 
builds on inconsistencies within the NSR approach in both 
the gap and the number equation. By far the most complete 
study of the NSR based theory for crossover was summarized 
in Ref. 8. By systematically introducing a series of improved 
approximations, the authors ultimately noted that one must in- 
corporate pairing fluctuation corrections into the gap as well 
as the number equation. 

It should be stressed that (with or without the approximate 
form for the number equation) the NSR scheme at T ^ 
was not designed to be consistent with the simple BCS-Leggett 
ground state, which they also discussed at length. This obser- 
vation was implicitly made elsewhere'*^ in the literature and 
can be verified by comparing the ground state density profiles 
based on the NSR scheme with those obtained in the Leggett 
mean field theory4^ It should also be stressed that T-matrix 
theories do not incorporate a direct pair-pair interaction; rather 
the pairs interact in an average or mean field sense. If one 
tries to extract the effective pairing interaction from any T- 
matrix theory, the absence of coupling to higher order Green's 
functions will lead to a simple factor of two relating the inter- 
boson and inter-fermion scattering lengths. More exact calcu- 
lations of this ratio lead to a factor of 0.6.''*'i2i^i^ 



C. Present T-matrix Scheme 

We now show that one obtains consistent answers between 
T-matrix based approaches and the BCS-Leggett ground state 
equations, provided the pair susceptibility contains one bare 
and one dressed Green's function. Thus, for simplicity, we 
refer to the present approach as "GGq theory". Throughout 



this paper we will emphasize the strengths of the present T- 
matrix scheme which rest primarily on a consistent treatment 
of pseudogap effects in the gap and number equations. This, 
in turn, leads to physical behavior for the thermodynamics, the 
superfluid density and the density profiles at all temperatures. 
Finally, we note that the present T-matrix scheme is readily 
related to a previously studied"**" approach to fluctuations in 
low dimensional, but conventional superconductors. A weak 
coupling limit of this GGq approach is equivalent to Hartree 
approximated Ginzburg-Landau theory.^'' 

We begin with the situation in which there is an equal spin 
mixture, and then generalize to the population imbalanced 
case. In the present formalism, for all T < Tc, the gap equa- 
tion is associated with a BEC condition which requires that 
the pair chemical potential ^ipair vanish. We will show below 
that because of this vanishing of Upair at and below Tc, to a 
good approximation one can move Go outside the summation 
in Eq. (01). As a result the self-energy is of the BCS-like form 



Thus 



<k 



(5) 



(6) 



Now we are in a position to calculate the pair suscepti- 
bility at general Q, based on Eq. (|2]l. After performing the 
Matsubara sum and analytically continuing to the real axis, 
ifti n + iO) we find the relatively simple form 



xiQ) = E [ 



/(i^k) - /(a-q) 

Ek - a-q + n + io+' 



(7) 



where u^,v^ — (1 ± Ck/-E'k)/2 are the usual coherence fac- 
tors, and f{x) is the Fermi distribution function. It follows 
that x(0) is given by 



X(0) 



2Ek 



(8) 



The vanishing of Upair (or generalized Thouless criterion) 
then implies that 

tpgi(0) = C/-i+x(0) = 0, T<Tc. (9) 

Substituting x(0) into the above BEC condition, we obtain the 
familiar gap equation 



2Ek 



(10) 



Here £'k = + A^, which contains the total excitation 
gap A instead of the order parameter Age- 

The coupling constant U can be replaced in favor of 
the dimensionless parameter, 1/kpa, via the relationship 
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TO/(47ra) = 1/U + X]k(2^fe)~^' where a is the two-body 
s-wave scattering length, and kp is the noninteracting Fermi 
wave vector for the same total number density. Therefore the 
gap equation can be rewritten as 



where 



m 
A'Ka 



E 



1 - V{Ek) 1 



2Eh 



(11) 



Here the "unitary scattering" limit corresponds to resonant 
scattering where a — > oo. For atomic Fermi gases, this scatter- 
ing length is tunable via a Feshbach resonance by application 
of a magnetic field and we say that we are on the BCS or BEC 
side of resonance, depending on whether the fields are higher 
or lower than the resonant field, or alternatively whether a is 
negative or positive, respectively. 

Finally, inserting the self energy of Eq. (|5]), into the Green's 
function, it follows that the number equation is given by 

n = 2Y}I{Ek)ul + f{-Ek)vl], (12) 



thus demonstrating that both the number and gap equation [see 
Eq. ( fT0l )1 are consistent with the ground state constraints in 
BCS-Leggett theory. 

Next we use this T-matrix scheme to derive Eq. (|5]l and 
separate the contribution from condensed and noncondensed 
pairs. The diagrammatic representation of our T-matrix 
scheme is shown in Fig. [T] The first line indicates the T- 
matrix, tpg, and the second the total self energy. The T- 
matrix can be effectively regarded as the propagator for non- 
condensed pairs. One can see throughout the combination of 
one dressed and one bare Green's function, as represented by 
the thick and thin lines. The self energy consists of two con- 
tributions from the noncondensed pairs or pseudogap {pg) and 
from the condensate (sc). There are, analogously, two contri- 
butions to the full T-matrix 



^ ^pg ^" ^sc : 

u 



(13) 



tsciQ) = ~^S{Q), (15) 

where we write Asc = —U X]k(^-kiCkt)- Similarly, we have 
for the fermion self energy 

E(i^) = + ^pg(K) = J2 t{Q)GoiQ - K) . (16) 

Q 

We see at once that 

^sc{K) = tsc{Q)GoiQ -K)^ -GohK)Al . (17) 
Q 

A vanishing chemical potential means that tpg{Q) diverges at 
Q = when T < Tc- Thus, we approximate*^ Eq. (fTSI l to 
yield 



A\T)^AliT) + A%iT), 
Importantly, we are led to identify the quantity Apg 



A 



iQ)- 



(19) 



(20) 



Note that in the normal state (where /ipair is non-zero), 
Eq. (fTsl l is no longer a good approximation. We now have a 
closed set of equations for addressing the ordered phase. We 
show later how to extend this approach to temperatures some- 
what above Tc, by self consistently including a non-zero pair 
chemical potential. This is a necessary step in addressing a 
trap as well.'*** 

The propagator for noncondensed pairs can now be quanti- 
fied, using the self consistently determined pair susceptibility. 
At small four-vector Q, we may expand the inverse of tpg, 
after analytical continuation, to obtain 



t-,\Q)^ain^ + z(n- 



2M* 



f-^pair 



(21) 



where below Tc the imaginary part Fg faster than as 
q ^ 0. Because we are interested in the moderate and strong 
coupling cases, where the contribution of the aifl^ term is 
small, we drop it in Eq. (|2T]) so that 



ipgiQ^ 



where we associate 



z- 



f-^pa 



Q 



2M* 



(22) 



(23) 



This establishes a quadratic pair dispersion and defines the 
effective pair mass, M*. This can be calculated via a small q 
expansion of x(Q)> 



1 



1 d'^x 



2M* 



Finally, one can rewrite Eq. ( |20] i as 

AlgiT)^z-'Y.f'ing), 



f2=0,q=0 



(24) 
(25) 



(18) 



where b{x) is the Bose distribution function. 

The superfluid transition temperature T^ is determined as 
the lowest temperature(s) in the normal state at which noncon- 
densed pairs exhaust the total weight of A-^ so that A^^ = A^ . 
Solving for the "transition temperature" in the absence of 
pseudogap effects'*'^-^"'^' leads to the quantity T^ '^ . More 
precisely, T^^^ should be thought of as the temperature at 
which the excitation gap A(T) vanishes. This provides a rea- 
sonable estimate for the pairing onset temperature T*. It is 
to be distinguished from T^, below which a stable superfluid 
phase exists. We note that T* represents a smooth crossover 
rather than a thermodynamic phase transition. 
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tpg \ + 



s = 



+ 



+ 



+ 




^sc 



Figure 1: T-matrix and self-energy diagrams for the present T- 
matrix scheme. The self-energy comes from contributions of both 
condensed (Esc) and noncondensed (Epg) pairs. Note that there is 
one dressed and full Green's function in the T-matrix. Here the T- 
matrix tpg can be regarded effectively as the propagator for the non- 
condensed pairs. 



It should be stressed that the dispersion relation for the non- 
condensed pairs is quadratic. While one will always find a lin- 
ear dispersion in the collective mode spectrum,^^ within the 
present class of BCS-BEC crossover theories, the restriction 
to a T-matrix scheme means that there is no feedback from 
the collective modes onto the pair excitation spectrum. In ef- 
fect, the T-matrix approximation does not incorporate pair- 
pair interactions at a level needed to arrive at this expected 
linear dispersion in the pair excitation spectrum. Nevertheless, 
this level of approximation is consistent with the underlying 
ground state wave function. 



The notion of phase separation between paired and un- 
paired states, separated by an interface, was first introduced 
in''° in the BCS limit, and it was more extensively discussed at 
T = in the crossover regime in Ref. 54 for the homogeneous 
case. A treatment of phase separation in a trap at zerc''^'''^ and 
at finite temperature lS lll6ll has received considerable recent at- 
tention. In a harmonic trap, phase separation leads to a nearly 
unpolarized gas at the center surrounded by a polarized, but 
essentially uncorrelated normal Fermi gas. Here one sees that 
the excitation gap A decreases abruptly to zero. By contrast, 
at higher temperatures, where the Sarma phase is stabilized, A 
decreases to zero continuously and there is a highly correlated 
mixed normal region separating a superfluid core and normal 
(uncorrelated) gas. 

We now extend the present GGq formalism to include po- 
larization effects. ''^ Including explicit spin indices, the pair 
susceptibility is given by 



1 



x(Q) = :^[xu(Q) + XiT(Q)] 



E 



i-KEk)-!{^,-k) 



f{Ek 



+ Ek — Cg-fe 



(26) 



where the coherence factors u^,v^ — (1 ± £,k/E\^) /2 are for- 
mally the same as for an equal spin mixture. For notational 
convenience we define 



III. GENERALIZATION TO INCLUDE POPULATION 
IMBALANCE 



It is relatively straightforward to include a difference in par- 
ticle number between the two spin species, within the con- 
text of the BCS-Leggett wave function. This is closely anal- 
ogous to solving for the spin susceptibility in BCS theory. 
The excitation energies are given by E)^^ = —h + Ek and 
Eki ^ h + Ek, where ^j. = ej, - ^ and E^ = ^/^^ + A^. 
Here /i — +/i|)/2 and h — {fi-^ — /X|)/2. We assume spin 
up fermions are the majority so that n-f > and h > 0. It is 
important to note that depending on h, fi, and A, the quantity 
i?fe| may on occasion assume negative values for a bounded 
range of fc-states. At T = this implies that there are regimes 
in fc-space in which no minority component is present. This 
leads to what is often referred to as a "gapless" phase. It was 
first studied by Sarma-^ at T = in the BCS regime. 

It is natural to extend this ground state Sarma or "breached 
pair" phase to include BCS-BEC crossover effects^^'^"*'^^'^^. 
The effects of finite temperatures were also studied using the 
current GGq, T-matrix scheme'***-^"-^^''''^, using the Nozieres 
Schmitt-Rink formalism^-^ as well as using an alternative 
many body approachi^i^ It should be noted, however, that 
the Sarma phase is generally not stable at T = except on 
the BEC side of resonance. Studies of the Sarma phase closer 
to unitarity and at low temperature reveal negative superfluid 
density^^ as well as other indications for instability.^** More 
generally, closer to unitarity, the Sarma phase stabilizes only 
at intermediate temperatures,^- while the ground state appears 
to exhibit phase separation. 



fix) = [f{x + h) + f{x~h)]/2, 



(27) 



Following the same analysis as for the unpolarized case, 
and using the above form for the pair susceptibility, the gap 
equation can be rewritten as 



Ana ^ 



l-2f{Ek) 
2Ek 



1 

2^ 



(28) 



The mean field number equations can be readily deduced 



Y.^f{Eko 



fi~Ek,)vl] 



where a = —a. The pseudogap equation is then 



A2^(T) = Z-i^6(f7,). 



(29) 



(30) 



Analytical expressions for Z and fig can be obtained via ex- 
pansion of x{Q) small Q (See, e.g., Ref. 62). This theory 
can readily be extended to include a (harmonic) trap as will 
be discussed in more detail in Sec. [Vll In case of a phase sep- 
aration, equilibrium requires T, /i, and and the pressure, p to 
be continuous across the interface or domain wall. Finally, it 
is useful to define polarization S in terms of 

N„{r) = j d\n,{r), N^N^^+N^, (31) 
5 = {N^-N{)IN. (32) 
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In this paper we do not discuss alternative phases such 
as the famous Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) 
states^ in which the condensate is associated with one or 
more non-zero momenta q. The competition between vari- 
ous polarized phases is associated''- with the detailed struc- 
ture of x{Q)- Indeed, there are strong similarities between 
these competing phases in polarized gases and Hartree-Fock 
theories which are used to establish whether ferro- or antifer- 
romagnetic order will arise in a many body system. The latter 
is associated with zero or finite wave-vector, respectively, and 
depends on the nature of the particle-hole spin susceptibility, 

^part^hole^Qy j^lis, in tUm, is given by ^P'^rt-hole^qs^ ^ 

+ Xo{Q), where Xo is the usual Lindhard function and U 
is the on-site repulsion. Here, by analogy the "ferromagnetic" 
case would correspond to the Sarma phase and the "antifer- 
romagnetic" situation to a LOFF like phase. Note, however 
that the relevant xiQ) necessarily involves the self consis- 
tently determined fermionic gap parameter A(T) and chem- 
ical potential /i, whereas for the magnetic analogue the bare 
particle-hole susceptibihty appears. 

IV. NORMAL-PHASE SELF-CONSISTENT EQUATIONS 

We next summarize the self consistent equations associated 
with the normal phase. We do not solve these at an exact level. 
This would require a numerical solution of the T matrix the- 
ory above Tc, which has been shown elsewhere^ to be very 
complicated. Instead we extend our more precise T < Tc 
equations in the simplest fashion above Tc, by continuing to 
parameterize the pseudogap contribution to the self energy in 
terms of an effective excitation gap A, using Eq. ( fTSl ), and 
thereby, ignoring the finite lifetime associated with the nor- 
mal state (pre-formed) pairs. We will, however make some 
accommodation of this lifetime in the following section. The 
self consistent gap equation is obtained from Eqs. ( |2TI ) and 
(O as 

t;g\0) ^ Zfipan- = U-' +xiO) (33) 

which yields 

^ + ^ = 2^Mpair , (34) 

k 

Similarly, above Tc, the pseudogap contribution to 

A2(r) = AL(r) + A2^(r) is given by 

q 

The density of particles can be written as 

n^2Y,KfiEk)+vlfi-Ek)] , (36) 

k 

It should be understood that the parameters appearing in the 
expansion of the T-matrix such as Z and ilq [See Eq. ( |22] |1 are 
all self consistently determined as in the superfluid state. 



In summary, when the temperature is above Tc, the order 
parameter is zero, and A — Apg. Since there is no conden- 
sate, fipair is nonzero, thus the gap equation is modified as 
= U^^ + x(0) ~ Z fipair- The number equation remains 
unchanged. From the above three equations, one can deter- 
mine /i, A and fipair- 



V. APPROXIMATE TREATMENT OF PAIR LIFETIME 
EFFECTS 

In the previous section, we discussed the extension of our 
more precise T < Tc equations above Tc, by continuing to 
parameterize the pseudogap contribution to the self energy in 
terms of an effective excitation gap A, using Eq. ( fTSl ), and 
thereby, ignoring the finite lifetime associated with the normal 
state (pre-formed) pairs. We will now make some accommo- 
dation of this lifetime by including "cut-off effects associ- 
ated with an upper limit of the momentum to be inserted into 
Eq. (O or Eq. (O. 

Below Tc, we can to a good approximation neglect the 
cutoff for the boson momentum q in evaluating the noncon- 
densed pair contributions to the pseudogap. This is justified 
by virtue of the divergence of tpg{Q) at Q — and low T 
so that the dominant contributions come from small q pairs. 
However, above Tc, pairs develop a finite chemical potential 
so that tpg{Q) no longer diverges and high momentum pairs 
would make substantial contributions to the integral in evalu- 
ating Apg via Eq. (|35] ). 

In order to make a more accurate evaluation, we take into 
account some aspects of the finite life time effects of the pairs. 
From Eq. O, one can read off the imaginary part as 

Imx(l^ + iO+,q) = Zro,q 

k 

+ [/(i?k) - f{^^-^)]vl5{E^ - a-q + f^), (37) 

where Fn q is the imaginary part of the pair dispersion. It is 
clear that Fn.q is nonzero when — min(i?k — Ck-q) < ^ < 
min(i?k + Ck-q) for any given (il, q). For on-shell pairs, we 
set = 57q — fjLpair in evaluating Ff^ q. Nevertheless, Ffj.q 
remains small for a large range of momentum q. Here we 
focus on positive pair dispersion so that the second term in 
Eq. ( |37] i vanishes. Apart from energy conservation imposed 
by the delta function, the factor 1 — /(i?k) ~ /(Ck-q) guar- 
antees that the contribution of the first term in Eq. dJTl l is very 
small when ^k-q < except at high T. As a very good es- 
timate, we impose a cutoff for q such that when q = qcut we 
hsNsilq — fipair = -Ek + ^k, wherc fc minimizes |^k|- Tokeep 
our calculations self-consistent, we also impose this momen- 
tum cutoff below Tc- 

At high enough T in the BCS and unitary regimes, we 
sometimes find that there is no solution for qcut when A be- 
comes small and — fipair becomes large. We then extrapolate 
qcut smoothly to zero at higher T via qcut oc \fK- This avoids 
the unphysical abrupt shut down of the pseudogap at high T. 
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Figure 2: 3D density profiles n(r) of a Fermi gas in a hiarmonic trap at unitarity at T/Tf ~ 0.01, 0.1, 0.2, and 0.3. The density distributions 
are smooth and monotonic, and become broader with T increasing. There is no bimodal feature in the density profiles, in agreement with 
experimental observations. Here Tf ~ Ef /ks is the global Fermi temperature and Rtf is the Thomas-Fermi radius. The density n(r) is in 
units of fc^"'. 



In the BEC regime, however, one finds that qcut = +oo and 
the pairs are bound and long Hved, as expected physically. 



VI. DENSITY PROFILES 

We now turn to include trap effects, with spherical trap po- 
tential Vexti'r) = \'muPr'^. Within a trap, we impose the 
force balance equation, — Vp = nWext, where p is the pres- 
sure and Vext is the trap potential. In the trap, the temperature 
is constant, so we have the relation Vp = nV/x. Thus we 
obtain V/x — —^Vext{r), or 

li{r)= no-Vext{r) (38) 

where = /i(0) and T4it(0) = 0. This shows that the force 
balance condition naturally leads to the usual local density ap- 
proximation (LDA) in which the fermionic chemical poten- 
tial ^ can be viewed as varying locally, but self consistently 
throughout the trap. 

We can readily extend our self consistent equations to in- 
corporate a trap, treated at the level of LDA. Tc is defined 
as the highest temperature at which the self-consistent equa- 
tions are satisfied precisely at the trap center At a temperature 
T lower than Tc the superfluid region extends to a finite ra- 
dius Rsc- The particles outside this radius are in a normal 
state, with or without a pseudogap. The important chemi- 
cal potential impair {r) is identically zero in the superfluid re- 
gion r < Rsc , and must be solved for self-consistently at 
larger radii. Our calculations proceed by numerically solv- 
ing the self-consistent equations. In the figures below, we 
express length in units of the Thomas-Fermi radius Rt f = 
y/2EF/imLu^) = 2(3iV)i/3/fcF; the density n{r) and total 
particle number N — J (fir n{r) are normalized by k\ and 
[kpRTpY ^ respectively. 

We determine Tc as follows: (i) An estimated initial value 
for chemical potential is assigned to the center of the trap 
/i(0), which determines the local /i(r) = /i(0) — Vcxtir)- (ii) 
We solve the gap equation ([U and pseudogap equation (l35l l at 
the center (setting Apg = A) to find Tc and A(0, Tc). (iii) We 
next determine the radius Rmax where A drops to zero, (iv) 



Next we solve the gap equation ( |34] | and pseudogap equation 
(|35T l for A(r, Tc) for r < Rmax- Then n{r) is determined us- 
ing Eq. ( [36] l. (v) We integrate n{r) over all space and enforce 
the total number constraint N — J (firn{r). We use nonlin- 
ear equation solvers which iteratively find the solution for the 
global /Lt(0) and the local gap parameters. Below Tc, an extra 
step is involved to determine the condensate edge, Rsc, where 
Asc drops to zero. Within the superfluid core, Eqs. ([T]l and 
(|35] | are solved locally for A and A^c, with pipair{r) = 0. 

A. Numerical results for unpolarized case 

In this section we address the particle density profiles at all 
T in the near-BEC, the near-BCS, and the unitary regimes. 
For the latter this work helps establish why the measured den- 
sity profiles appear to be so featureless. ^^'^^ Some time ago 
it was found^^ that at unitarity the profiles were reasonably 
well described by a Thomas-Fermi (TF) fit at zero T, and in 
recent work^^ this procedure has been extended to finite tem- 
peratures, suggesting that it might be quite general. Our cal- 
culations indicate this TF fit is reasonably good below Tc, and 
becomes substantially better above Tc- The width of the pro- 
files has been used to extract an effective temperature scale.— 
If we follow the same procedure^** on our theoretical profiles 
we find that the temperature scale coincides with the physical 
T quite precisely above Tc- Below Tc, because the conden- 
sate edge moves inwards as temperature increases, this tends 
to compensate for thermal broadening effects. In this way, 
in the superfluid phase the effective temperature needs to be 
recalibrated^^ to arrive at the physical temperature scale. 

Our work differs from previous theoretical studies "is^ by 
including the important effects of noncondensed pairs^i^l 
which are associated with pseudogap effects. These "bosons" 
are principally in the condensate region of the trap, whereas 
fermionic excitations tend to appear at the edge where the gap 
is small. In contrast to the work of Refs. |42| andfsl, our den- 
sity profiles are monotonic in temperature and show none of 
the sharp features in the BEC which were predicted^** from a 
generalization of the Nozieres-Schmitt-Rink approach. Our 
calculations show that pseudogap effects are responsible, not 
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Figure 3: Comparison of 3D density profiles n{r) at different temperatures between the unitary (left), near BEC (l/kpa = 1, middle) and 
near BCS (l/kpa = —0.5, right panel) regimes. They broaden with increasing T but shrink with increasing pairing strength. 



only for the relatively featureless density profiles we find in 
the unitary regime, but also for the behavior of the associated 
temperature evolution. 

Figure |2] shows the behavior of the three-dimensional (3D) 
density profiles of a Fermi gas at unitarity as temperature pro- 
gressively increases (from left to right). One can see that the 
profiles become progressively broader with increasing T. Be- 
cause there is no bi-modality or other reflections of the con- 
densate edge, one can thereby understand why the Thomas- 
Fermi fits are not inappropriate. A more quantitative compar- 
ison of this unitary case with experiment is in Ref. 70. 

In Fig. Owe present a comparison of the density profiles in 
a unitary system with the near BEC and near BCS cases. On 
the BEC side of resonance (l/kpa = 1) the profile is signif- 
icantly narrower than that on the BCS side. The unitary case 
is somewhere in between. The quantity (3 which is used in the 
literature to parameterize this width is of the order of —0.41 
as compared with experiment where f3 « —0.55. Conven- 
tionally, f3 is defined as the ratio of the attractive interaction 
energy to the kinetic energy and is given hy fi — {1 + P)Ep 
and /io = \/l + /3Ef for homogeneous and trapped unitary 
gases, respectively. The discrepancy between theory and ex- 
periment is associated with the absence of Hartree self-energy 
corrections in the BCS-Leggett mean field state. Thus, for 
more quantitative comparison with unitary experiments^ we 
match the /? factor by going slightly on the BEC side of reso- 
nance. 



B. Numerical results for polarized case 

In this section we show how the general shape of the den- 
sity profiles at unitarity changes as one varies the polarization. 
Unlike the unpolarized case, we can identify features in the 
polarized gas profiles which indicate whether or not the gas 
is superfluid; these features are rather similar to what is ob- 
served experimentally .^ i,7 1 .72.73 ^j^^ trace the evolution of 
the profiles from phase separation at low temperature to the 
Sarma phase. 

We begin with Fig. |4] which shows the phase diagram at 
l/kpa = 1.5 on the BEC side of resonance. This should be 
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Figure 4: Phase diagram in the T-S plane of a population imbalanced 
Fermi gas at l/kpa — 1.5 in the BEC regime. Here PS denotes 
phase separation, which exists only at low T and high polarization. 



compared with the counterpart phase diagrams for unitarity 
and the near-BCS which have been presented in Ref. 61. The 
principal difference between unitarity and this case is that for 
the former the phase separation (PS) region is present at low 
T over the entire range of polarizations, whereas in the BEC 
regime, it has been pushed toward the high polarization region 
of the phase diagram. 

Figure |5] shows the temperature dependence of the unitary 
profiles for majority and minority spin components at p = 0.5 
for a range of temperatures, increasing from left to right. The 
lowest temperature {T/Tp = 0.01) corresponds to a situa- 
tion when phase separation is present, while the three higher 
temperature correspond to the Sarma phase. The condensate 
edge is clearly apparent in the phase separation scenario, with 
a jump in order parameter at the edge. For the Sarma phase 
cases, bimodality is clearly visible in the minority profile, and 
a kink-like feature is present in the majority profile well be- 
low Tc- At high T, both majority and minority profiles be- 
come closer to a Thomas-Fermi distribution, as polarization 
has penetrated into the superfluid core. 

The vertical dashed lines for the three Sarma cases in the 
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Figure 5: Evolution of the 3D density profiles n(r) with temperature T at unitarity and polarization S — 0.5. The upper (black) and lower (red) 
curves are for the majority and minority species, respectively. From left to right, T/Tf = 0.01, 0.1, 0.2, and 0.3. Phase separation occurs for 
T/Tf ~ 0.01, where the profile shows abrupt changes at the phase boundary, whereas the Sarma phase prevails in the other cases so that the 
profiles are smooth. Bimodality is clearly visible in the Sarma cases. Within the vertical dashed lines there exists a paired normal region. 



figure delimit the paired normal region. They correspond to 
the condensate edge, where A^c drops to 0, and the gap edge 
where the total excitation gap A = Ap^ smoothly disappears. 
Between the two dashed lines the system is in a paired or 
highly correlated mixed normal state. ^''^^ The width of this 
mixed normal region grows with increasing temperature, and 
the condensate edge disappears above Tc- Outside the gap 
edge, the gas is free; there is a small range of r where both 
spin components are present and a wider range where only 
the majority appears. In the phase separation regime, such a 
mixed normal region is essentially absent^^ii^^ and the con- 
densate edge is indicated by a single dashed vertical line. For 
low T, we note that the condensate is essentially unpolarized. 

In summary, in the phase separation regime, there are sharp 
discontinuities in the profile associated with the condensate 
edge, the other side of which is a free Fermi gas. In the Sarma 
phase, which is stabilized at higher T, there may also be in- 
dications of the condensate edge. Beyond the superfluid core, 
there is a highly correlated mixed normal region which caiTies 
a significant fraction of the polarization and is associated with 
the pseudogap phase. Finally, in the outer regime of the pro- 
file there is a free Fermi gas, which may consist of majority 
only or of both spin states. These three regions in the Sarma 
phase seem to be in accord with experiment.^'-^'' An impor- 
tant additional finding is that except at high temperatures the 
superfluid core seems to be robustly maintained at nearly zero 
polarization, as observed experimentally.A'iZii22iS 



VII. THERMODYNAMICS 

In this section we introduce^"^ an approximate form for the 
thermodynamical potential (density), fl. We can, to a high 
level of accuracy, write this down analytically. It is impor- 
tant to assess this approximate form by studying various ther- 
modynamical identities. We will do so here by checking 
Maxwell's relations as well as establishing the relationship 
p — between energy density E and pressure p, which 
is expectediSiii to apply at strict unitarity. In the superfluid 
phase, we find there is essentially no deviation from the pre- 
cise thermodynamical relations. Above Tc, we find deviations 



of from one to a few percent. 

We begin with the unpolarized case. The quantity fl is asso- 
ciated with a contribution from gapped fermionic excitations 
ilf as well as from non-condensed pairs, called flf,. These 
two contributions are fully inter-dependent. The gap in the 
fermionic excitation spectrum is present only because there 
are pairs and conversely. We have 

nf = A2xo + 5][(a-£;tc)-2rin(l + e-^''/^)], 



where xo 



q 



(39) 



Zp-pair- The pressure is simply 

p^-n (40) 



Here ppair = at T < Tc, while above Tc the superconduct- 
ing order parameter Age — 0. Providing that we ignore the 
very weak dependence of the parameter Z and the pair mass 
M* on A, /i and h, we are able to derive our self consistent 
gap, pseudogap and number equations variationally. These 
self-consistent (local) equations are given by 



on 

dA 







(41) 



which represents the gap equation ( |34| |. Similarly, we have 







(42) 



which leads to the equation for the pseudogap given by 
Eq. ( [35l ). Finally, the number equation 



(43) 



which yields Eq. Q6j. In a trap, this is subject to the total 
number constraint = J (Prn{r). 
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From the above thermodynamical potential, we can deter- 
mine all other thermodynamic quantities. The energy (den- 
sity) is 

E = Ef+Eb 
Ef = -^\o + Y.'^{i^.-E^)-2EU{E^.)]+|J^n, 



Eb = ^ f7q6(f7q — flpair) 

q 

and the entropy (density) is 

S = Sf + Sb 



(44) 



ln(l 



Sb 



E 



^6(r!q)+ln(l-e-^^'' 



(45) 



It is easy to verify the relation 

VLf ^ Ef -TSf - (46) 

and 

^b Eb T Sb f^pair'^pair (47) 
with Upair = ZA^g. 

In the actual calculations of thermodynamic properties we 
combine Eq. ( |39] l with a microscopic calculation of the non- 
condensed pair propagator, thereby determining Z, and ilq 
from the expansion of the inverse T-matrix. We test the valid- 
ity, then, of our expression for the thermodynamic potential 
51 by examining Maxwell identities. Indeed the deviation is 
generally at most at the few percent level, as will be illustrated 
below. 

Finally, we end our analytical discussion with expressions 
for a polarized gas. Here the thermodynamical potential is 
given by 



n = 0/ + r^b 



^rin(i 



nb 



^Tln(l 



(48) 



Competing with this phase is the free Fermi gas phase which 
has thermodynamical potential density 



free 



-r^ln(l + e-« 



r/T 



(49) 



Here Ey^^ = E]^ ^ h and ^kcr = £.k T h for spin a =t, i, 
respectively. 

It should be noted that in this paper, we are concerned with 
primarily the internal energy (density) and pressure without 
the contribution from the external trap potential, in order to 
test the relationship p/ E — 2/3. The internal energy can be 
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Figure 6: Thermodynamic behavior of a homogeneous Fermi gas at 
different l/kpa as labeled. Shown in (a) the comparison between 
per-particle energy (multiplied by 2/3, dashed lines) {2/3)E/N and 
pressure p/n (solid lines) and in (b) the entropy per particle S/Nks- 
Here N = n since we have set volume V = 1. 



obtained by substituting for the chemical potential the local 
/i(7-) in the term Ef in Eq. (l44l i. The total energy, which in- 
cludes the trap potential, may be obtained by further adding 
nVext{f) to Ef in Eq. ( l44b . For a harmonic trap at unitarity, 
the internal energy and the external trap potential energy are 
equal.— 



A. Numerical Results for unpolarized case 

In this section we discuss numerical results for thermody- 
namic properties principally for trapped Fermi gases within 
the unitary, near-BCS and near-BEC regimes. We find that 
unpaired fermions at the edge of the trap, where A is small, 
provide the dominant contribution to thermodynamical vari- 
ables such as E and S at all but the lowest T. In addition 
to the usual gapped fermionic excitations, there are "bosons" 
which correspond to finite momentum pairs. Above Tc these 
"bosons" lead to a normal state fermionic excitation gap (or 
"pseudogap")i2i2iI^i2^ They are dominant only at very low 
T <C Tc, leading to 5 oc T^/^. We emphasize that the normal 
state of these superfluids is never an ideal Fermi gas, except 
in the extreme BCS limit, or at sufficiently high T above the 
pseudogap onset temperature T*. 

In Fig.|6l we plot (a) the energy per-particle E/N (dashed 
lines) multiplied by 2/3 and pressure p/n (solid lines) and (b) 
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Figure 7: Test of Maxwell relations. The solid and dashed curves are 
dn/dT and ds/d/i, respectively, as function of trap radius, at dif- 
ferent temperatures for l/kpa = (upper) and l/kpa = 1 (lower 
panel). As labeled, the black, red and green colors correspond to 
T/Tf ~ 0.01, 0.15, and 0.3, respectively. The difference between 
the solid and dashed curves, while largest in the normal regime, is 
almost negligible. 
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Figure 8: Trap averaged (per-particle) pressure pfN (solid), and in- 
ternal energy Ej^ in the upper panel and entropy SjN in the lower 
panel as a function of T for l/fcFC[ = — 0.5, 0, and 1, as labeled. The 
relation p = 2E/?> is satisfied for and only for the unitary case. The 
agreement is nearly perfect at T < Tc- In the pseudogap phase, 
the discrepancy remains very small (< 5%). Here the contribu- 
tion from the external trap potential is not included in the p or E. 
Tc/Tf = 0.19, 0.28, and 0.33, respectively, for the three regimes. 



the entropy S/Nkg for a homogeneous system and for a range 
of values of l/kpa, from noninteracting (1/kpa = — oo) to 
near BEC {l/kpa = 1/2). It can be seen that all curves ap- 
proach the free Fermi gas results at T > T*. It is also clear 
that, as expected, the energy and entropy are lowered as the 
system goes deeper into the BEC. The pairing onset temper- 
ature T* stands out in the figure as the most apparent tem- 
perature scale. We find virtually no thermodynamic feature 
at Tc- A small feature should be present in the BEC, becom- 
ing larger as the BCS regime is approached. This would ap- 
pear if we included lifetime effects associated with the non- 
condensed pairs; in order to make the calculations manage- 
able, we have ignored this complexity which has been ad- 
dressed elsewhere. It should be stressed that T* represents a 
crossover temperature and is not to be associated with singular 
structure in thermodynamical variables, unlike Tc- 

The comparison between the dashed and solid lines in 
Fig. |6ja) represents an important indicator of the universal- 
ity expected at strict unitarity, where the energy density and 
pressure satisfy p — Indeed the two curves are virtu- 

ally indistinguishable in the superfluid phase at unitarity, and 
remain very close to each other in the normal phase. This rela- 



tionship also holds for the non-interacting gas. By contrast, on 
the BEC side of resonance this relation is seriously violated, 
as expected. 

Figure [T] represents a test of one particular Maxwell rela- 
tion for the unitary case (upper panel) and for the near-BEC 
(l/kpa — 1, lower panel). Here we compare dn/dT (solid 
lines) with ds/dfi (dashed lines). The horizontal axis is the 
trap radius in units of Rtp- At the lowest temperature this 
Maxwell relation is very well satisfied. The feature shown 
in the plotted derivatives corresponds to the condensate edge. 
As the temperature is raised the deviation is slight, but per- 
ceptible. The small breakdown in the Maxwell relations cor- 
responds to our approximate treatment of the normal phase as 
discussed in Sec. lIVI 

In Fig. [8] we plot the trap averaged pressure (per particle) 
p/N (solid) and {2/3)E/N (dashed lines) in the upper panel 
as well as entropy S/N in the lower panel, as a function of 
temperature. For each quantity, the three curves correspond 
to unitarity and near-BCS (1/kpa = —0.5) and near-BEC 
(1/kpa = 1), respectively, as labeled. As for the homoge- 
neous case in Fig.|6] the closer the system is to BEC the lower 
the energy and entropy, as expected. Although not shown 
here, all curves will approach the free Fermi gas curve at suffi- 
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Figure 9: Spatial profile of tlie pressure p{r) (solid) and internal en- 
ergy (2/3)i?(r) (dashed lines) in the main figures as well as the en- 
tropy S{r) in the insets at T /Tf = 0.01 (black), 0.15 (red), and 0.3 
(green curves), for the unitary (upper panel) and near BEC (lower 
panel), respectively. The arrows point in the direction of increasing 
T. At unitarity, the relation p = 2E /3 is nearly perfect for the low 
T profiles, while in the pseudogap phase, the deviation is less than 
5%. The l/kpa — 0.5 case clearly violates the p — 2E/3 relation. 
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Figure 10: Temperature dependence of the trap averaged pressure 
p/N (solid), internal energy E/N (dashed) in the upper panel and 
entropy S/N in the lower panel in a trap at unitarity for polariza- 
tion (5 = 0.1 (black), 0.5 (red), and 0.8 (green curves), as labeled. 
The p — 2E/3 relation is satisfied at unitarity even with population 
imbalance. The small kink in S/N indicates the change from phase 
separation to Sarma state. Tc/Tf = 0.28, 0.25, 0.19 for polarization 
5 = 0.1, 0.5, and 0.8, respectively. 



constant dependence on the trap radius. 



ciently high T, corresponding to their respective T*. By com- 
paring the solid and dashed lines in the upper panel, one can 
see that the relation p = 2E/3 is essentially satisfied at uni- 
tarity. 

Figure |9] plots the spatial distribution of the pressure p 
(solid) and the energy 2E/3 (dashed lines), as well as the en- 
tropy S (inset) for three different temperatures, for the unitary 
(upper panel) and the near BEC (l/fci?a = 1/2, lower panel) 
cases, respectively. The relation p/E — 2/3 holds very well 
at unitarity for all temperatures shown, but, as expected, it 
is clearly violated in the near BEC case. For 1/kpa = 1/2, 
one sees that the energy becomes negative at intermediate trap 
radii. This reflects the fact that at these radii, the density is re- 
duced so that the local quantity 1 /kpa is effectively increased 
and the gas is in the BEC regime. At unitarity the entropy 
in the inset tends to peak towards the trap edge; this reflects 
the contribution from free fermions. By contrast these free 
fermions are relatively absent in the near-BEC case and the 
entropy is dominated by pair excitations leading to a relatively 



B. Numerical results for polarized case 

In this section we discuss the behavior of thermodynam- 
ical variables for a polarized gas at unitarity. In the upper 
panel of Fig. [TO] we compare the trap averaged pressure per 
particle, p/N (solid curves) and energy {2/3)E/N (dashed 
curves) as a function of temperature, for three different po- 
larizations 5 — 0.1, 0.5, and 0.8. The lower panel shows the 
corresponding behavior of the entropy S/N. The figure illus- 
trates that the lower the polarization the lower is the energy 
and entropy. This is because the system can take full advan- 
tage of the pairing when the polarization is small. Importantly, 
the upper panel demonstrates that the relation p/E = 2/3 also 
appears to hold for a polarized gas. There are small kinks in 
the entropy curves at the two higher polarizations which re- 
flect the transition from the phase separated to Sarma state. 

The spatial profiles of the three thermodynamical variables 
are plotted for three different temperatures in Fig.[TT]at fixed 
polarization 6 — 0.5. The results are not dramatically dif- 
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Figure 11: Spatial profiles of pressure (solid) and energy (dashed) 
in the main figure and entropy in the inset at unitarity and polariza- 
tion 5 = 0.5, for T/Tf = 0.01 (black), 0.15 (red), and 0.3 (green 
curves), respectively. The black, red and green curves correspond to 
T = 0.01, 0.15, 0.37f, respectively. The arrow points along the di- 
rection of increasing T. The p — 2E/3 relation is essentially perfect 
for low T profiles, and the deviation remains very small (< 5%) in 
the pseudogap phase. 



ferent from the unpolarized case shown in the upper panel of 
Fig lH One can see that the p/E = 2/3 relation holds rather 
well across the trap and that at intermediate temperatures, the 
entropy tends to peak somewhat inside the trap edge, reflect- 
ing the excitations of nearly free fermions in this regime. 



VIII. SUPERTLUID DENSITY 

An essential component of any theory for BCS-BEC 
crossover is establishing that the superfluid density is well be- 
haved. The superfluid density (T) is perhaps the best re- 
flection of a proper (or improper) description of the superfluid 
phase. This meaningful description is not at all straightfor- 
ward to come by once one includes self energy corrections to 
the BCS gap and number equation. These two must be treated 
on an equal footing in order for the "diamagnetic" and "para- 
magnetic currents" to precisely cancel at Tc when approached 
from below. (And the Tc that one computes from below has 
to be the same as that computed from the pairing instability of 
the normal phase). 

This cancellation of diamagnetic and paramagnetic currents 
is deeply and importantly related to generalized Ward iden- 
tities as we will show below. These arise from a connection 
between the one particle properties (which show up in the dia- 
magnetic current, through the number equation) and the two 
particle properties (which, for example, reflect the fermionic 
excitation spectrum and show up in the gap equation). It 
is important to stress at the outset that because we must dis- 
tinguish between the gap and the order parameter, there is no 
unambiguous way to make use of the Nambu Gor'kov formal- 
ism. One can readily see, however, that the combination GGq 



is, in effect, proportional to that Gor'kov "F" function which 
involves the full excitation gap A, rather than the order pa- 
rameter 

Whether one considers a charged or an uncharged system, 
the formal analysis is the same. Here for the sake of defi- 
niteness we refer to a charged superconductor. We consider 
the in-plane penetration depth kernel K{0) in linear response 
theory. Within the transverse gauge we may write down this 
response without including the contribution from collective 
modes. The London penetration depth is A^^ = ^oe^{ns/m), 
where /io is the magnetic permitivity. Here we set /ig = e = 1 
for convenience. From linear response theory. 



(0) - (-) - PxxiO) , (50) 
\m/ XX 



where K is defined by 



(51) 

and the current-current correlation function 



l3 



= / e*"""OM(q, ^)>(-q> 0)) 



(52) 



= -2 ^ A^{K, K + Q)G{K + Q)K{K + Q, K)G{K) . 

K 

Here we use the four-vector notation, = (0, A), = 
and the bare vertex = (1,A). Summation is as- 
sumed on repeated indices, with the convention A^B^ = 
A^Bo — A • B. Without loss of generality we can ignore 
collective mode effects and work in a transverse gauge. 
For the bare vertex, we have Aq = 1 and 



X{K, K + Q)^ Vkek+q/2 = - (k + § 
' m V 2 



(53) 



The electromagnetic vertex can be written in terms of the cor- 
rections coming from the two self-energy components as 



A = A + (5A.„„ + 5Ks 



(54) 



where 5Kpg is the pseudogap term. This contribution deriving 
from pair fluctuations contains terms associated with Maki- 
Thompson (MT) like diagrams as well as Aslamazov-Larkin 
terms (AL) which appear in the theory of conventional super- 
conducting fluctuations. Here the situation is somewhat more 
complex because of the appearance of one dressed and one 
bare Green's function in the pair propagator, which leads to 
two AL diagrams. As a result the AL term itself depends on a 
(gauge covariant) vertex function A'. We may write 



SKpg EE SKmt + Mali + SKalA^') 



(55) 



The diagrams contributing to the full electromagnetic vertex 
A in the transverse gauge are given in Fig. [12] Here Amt is 
given by the MTpg diagram, and SAgc is given by the MTgc 
diagram. In contrast to the electromagnetic vertex A, the 
gauge covariant vertex A' satisfies a generalized Ward iden- 
tity to be discussed below. 
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Figure 12: Diagrams contributing to the full electromagnetic vertex 
A in the transverse gauge. Here the wiggly lines represent the T- 
matrix tpg and the dashed line represent the singular "condensate 
propagator" tsc, both shown in Fig.[T] The gauge covariant full ver- 
tex A' contains the electromagnetic vertex insertion along tsc- 



We now show that there is a precise cancellation between 
the MTpg and ALi pseudogap diagrams at Q = 0. This 
cancellation follows directly from a generalized Ward iden- 
tity (GWI) 

Q ■ X{K, K + Q) = G^\K) ~ Gq\K + Q) , (56) 
which can be shown to imply 

Q ■ [SAal, [K, K + Q)+ SAmt{K, K + Q)] = Q (57) 

so that SAali {K, K) = —SAmt{K, K) is obtained exactly 
from the Q — > limit of the GWI. 
To see this explicitly note that 



-J2tiP)Go{-K-Q + P) 
p 

xA^(-if -Q + P,-K + P) 
xGo{-K + P) (58) 



Similarly we have 

^Kl, - -J2Goi-K + P)tiP + Q) 
p 

X { ^ Gi-K' + P)Go(A" + Q) 

K' 

xX^iK' + Q,K')Go{K')}t(P) (59) 

We may write 



tiP)-' = U-' -J2GiKi+P)Goi-Ki) (60) 
Then combining terms 



MT 



SA 



ALi 



t{P)Goi-K + P) 



X [Go{-K -Q + p)[g^Hp -K)- G^\P -K-Q)] 
- t{P + 0) ^ G{-K' + P)Go(X' + Q)Go(X') 

K' 

x[Gi,\K')-Gi,\K' + Q)]] 

It then follows using Eq. ( l60l ) that this equation vanishes 
and we have proved the desired relation between the Maki- 
Thompson vertex and the ALi vertex. 

The GWI is not to be imposed on A since we are evaluat- 
ing the electrodynamic response in a fixed (transverse) gauge. 
However, the full gauge covariant internal vertex A' is consis- 
tent with the GWI. This internal vertex A' then satisfies 

Q-A'{K,K + Q) = G-\K) - G-\K + Q) . (61) 

The above result can be used to infer a relation analogous 
to Eq. (|52| for the AL2 diagram; so that 6Aal2 (K, K) = 
SAmt{K, K). More generally 

Q ■ {5Aali + SAal2) = -2Q ■ SAmt, (62) 

Therefore the combination of these three diagrams (in con- 
junction with Eq. ( fSSl l) leads to 

Q ■ 6Apg{K, K) = -Q- SAmHK, K) , (63) 

which expresses this pseudogap contribution to the vertex en- 
tirely in terms of the Maki-Thompson diagram shown in the 
figure. One can show explicitly that 



5A^l„{K,K) = - 



-'pg 
dkn 



(64) 



This can be proved as follows. We write 

Q-SAmt = -Y,tp9iP)[Go{-K+P)-Go{-K-Q+P)], 
p 

(65) 

where we have used the GWI involving the bare Green's func- 
tions to eliminate A. Now taking the q = limit with oj — 
and using Eq. ( |63] l and the expression of Epg(A') we arrive at 
Eq. dSl- 

Combining terms we find 



5A';;JK,K) = 



(66) 



This demonstrates consistency; that is, the usual Ward identity 
applies to the pseudogap contribution. 

Now we turn to the superconducting vertex contributions. 
As can be seen by a simple inspection of the diagrams, the 
superconducting contribution is closely analogous to Eq. (|64] i 
so that we have 



16 




Figure 13: Spatial profiles of superfluid density at zero polarization at different temperatures (as labeled) for l/kpa — (left), 1 (middle), 
and -0.5 (right panel). The insets show the T-dependence of trap integrated superfluid density. All the profiles are smooth, single- valued, and 
monotonic, evolving continuously with radius and temperature. 



More generally, we can define a relationship 



(67) 



Importantly, the above equation contains a sign change (as 
compared with Eq. (l66ll). This is associated with the trans- 
verse gauge and violates the Ward identity. It is central to the 
existence of a Meissner effect. The fact that the pseudogap 
contributions are consistent with generalized Ward identities 
is an important aspect of the present calculations. This im- 
plies that there is no direct Meissner contribution associated 
with the pseudogap self-energy. 

We next explicitly evaluate the superfluid density using 
Eq. dSOb . For this purpose, we only need the spatial compo- 
nents of the vertex functions. Note that the pseudogap contri- 
bution to [ns/m) drops out by virtue of Eq. (l66l ). The density 
can be rewritten using integration by parts. 



dek dG(K) 



dka 



K 



dka \dkf3 dk/s 



dkp 

' dkp J ' 
(68) 



where a, (3 — 1,2,3. Note here the surface term vanishes in 
all cases. The superfluid density is given by 



K ^ 



SAsc{K, K),- 



dkx 



(69) 



Equation ( |69l ) can be readily evaluated using the su- 
perconducting vertex and the superconducting self-energy 
Tjsc{K) — — A^^Go(— A') associated with our GGo-based T- 
matrix approach. In addition, we introduce an approximation 
in our evaluation of G via Eq. ( fTsT i, to find 



2 

^ El 



1 - 2/(Sk) 



2Eu 



\dkj 



(70) 



\m J \m J 



BCS 



(71) 



where (ns/m)^^ is just {ns/m) with the overall prefactor 
A^^ replaced with A'^ in Eq. dTOl i. Obviously, in the pseudo- 
gap phase, {ris/m)^'^'^ does not vanish at T^. 

Finally, in the polarized case it can be shown that the super- 
fluid density is given by Eq. ( iTOl i with the Fermi function and 
its derivative replaced by the quantities / and /', respectively. 



A. Numerical results for unpolarized and polarized cases 

The behavior of the superfluid density Ug (T) is viewed as 
one of the important indicators of the quality of a given BCS- 
BEC crossover theory. Plots of ns{T) in Ref. 78 stop at about 
Tc/2, above which it is argued that the calculations are un- 
reliable. Alternative plots^** show double-valued functions, 
particularly on the EEC side of resonance. While ns{T) is 
not explicitly evaluated, it will necessarily exhibit a first order 
transition in the work of Ref. 

It is important, then, to show that (T) corresponds to the 
appropriate physical behavior in the current theory. First, we 
present results for unpolarized Fermi gases. The spatial dis- 
tributions of ns{r) in a trap are plotted in Fig. [13] for differ- 
ent temperatures and three different scattering lengths ranging 
from near BCS to unitary to near EEC. In the insets are plot- 
ted the temperature dependence of the trap integrated super- 
fluid density. All curves are well behaved, single-valued, and 
monotonic from T — to T — Tc- The superfluid density 
vanishes precisely at Tc- 

Analogous plots are shown in Fig. [14] for a polarized gas 
in the unitary case and at three different polarizations S = 
0.1, 0.5, and 0.8. The main figures present plots as a function 
of trap radius, whereas the insets are plots as a function of 
temperature. Here, by contrast, the behavior is not always 
smooth. These sharp features are all expected and associated 
with polarization effects. At the lowest temperatures in the 
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Figure 14: Spatial profiles of superfluid density at unitary at different temperatures (as labeled) for polarization 5 = 0.1, 0.5, and 0.8 from 
left to right. The insets show the trap integrated superfluid density as function of T. High T profiles are in the Sarma phase, and therefore, 
smooth, evolving continuously with radius. In contrast. The lowest T curves are in the phase separation regime and thus show an abrupt drop. 
The kinks in the trap integrated Us reflect the transition from phase separation to Sarma state. 



main body of each of these figures one can see the effects of 
phase separation on Ug. The superfluid density stops abruptly 
at the interface between the normal and superfluid. At higher 
T in the Sarma phase, the curves end continuously at the trap 
edge. At the higher two polarizations the two insets indicate 
kinks which reflect the transition from a phase separated to a 
Sarma phase. 



IX. CONCLUSIONS 

There are many different renditions of BCS-BEC crossover 
physics in the literature, but what has guided us here is the 
implementation of a sound methodology for characterizing 
three fundamental properties: thermodynamics, density pro- 
files and superfluid density with and without population im- 
balance. While there is considerable emphasis in the litera- 
ture on numerical precision one goal of this paper was to set 
up a different set of criteria against which theories as well 
as simulations can be checked. Monte Carlo simulations are 
sometimes argued**" to be the ultimate theory. While they may 
provide reliable numbers, these alone (in the absence of more 
analytic many body schemes) will not yield sufficient insight 
into the complex physics of these very anomalous superfluids. 

Four important and inter-related physical properties were 
emphasized here, (i) There must be a consistent treatment of 
"pseudogap" effects. That is, the fermionic excitation spec- 
trum, ii^k must necessarily be modified from the usual BCS 
form. Here, based on a systematic analysis, we implement 
this modification by replacing the order parameter with the 
total excitation gap A. (ii) The theory must lend itself to a 
consistent description of the superfluid density ns{T) from 
zero to Tc. The quantity ns{T) should be single valued and 
monotonic.^ It must necessarily disappear at the same Tc one 
computes from the normal state instability; ns{T) is at the 
heart of a proper description of the superfluid phase, (iii) The 
behavior of the density profiles, which are at the basis for all 
thermodynamical calculations of trapped Fermi gases, must 
be compatible with experiment. Near and at unitarity, they are 



relatively smooth and featureless, well fit to a Thomas-Fermi 
like form. Only in the presence of polarization effects can 
one use these unitary profiles to find signatures of the conden- 
sate edge, (iv) The thermodynamical potential Q. should be 
variationally consistent with the gap and number equations. It 
should satisfy appropriate Maxwell relations and at unitarity 
be compatible with the constraint relating the pressure p to 
the energy density: p = Here we find this to be the case 
for a population imbalanced gas as well to the same level of 
numerical precision as for an unpolarized gas. 

For semi-quantitative comparisons with experiment there 
have been notable successes within the present theoretical 
framework which address a very wide group of experiments, 
including polarized and unpolarized gases. '*'^'-^^'™''*' **^'^^ 
However, it is clear that detailed quantitative agreement is 
not always possible. The calculated (3 factor at unitarity 
{j3 = —0.41), is not precise, as compared with experiment 
{f3 w —0.55). Moreover, the ratio of effective inter-boson 
scattering length to the fermionic scattering length is found 
to be 2.0, rather than 0.64^ Indeed, inter-boson effects are 
included only in a mean field sense at the level of the sim- 
ple BCS-Leggett wave function and related T-matrix scheme. 
One knowsii how to arrive at a more Bogoliubov-Uke treat- 
ment of the pairs which properly treats inter-boson effects 
appropriate to the deep BEC. It can be shown'"* to yield the 
factor 0.6. This involves adding to the wave function addi- 
tional terms involving four and six creation operators. How- 
ever, there is no natural and tractable extension at unitarity. 

We have emphasized here that what is most unique and in- 
teresting about these trapped Fermi gases lies not so much in 
the ground state, but rather in finite temperature phenomena. 
It is at finite T that one sees a new form of fermionic superflu- 
idity in which pair condensation and pair formation take place 
on distinctly different temperature scales. This temperature 
separation requires radical changes in the way we think about 
fermionic superfluidity, relative to our experience with strict 
BCS theory. We have argued here that at this relatively early 
stage of our understanding, it is more important to capture the 
central physics of this exotic superfluidity, than to arrive at 
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precise numerical agreement with experiment. Ultimately we 
must do both, as has been possible for the Bose gases. Never- 
theless assessing a theory based on understanding the qualita- 
tive physics has to proceed an assessment based on quantita- 



tive comparisons. 
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